from math import *

TWOPI=pi*2.0
PIOVER2 = pi/2.0
DE2RA = 0.01745329252
RA2DE = 57.2957795129
ERAD = 6378.135
ERADM = 6378135.0
AVG_ERAD = 6371.0
EPS = 0.000000000005
FLATTENING = 1.000000 / 298.257223563 # Earth flattening (WGS84)
KM2MI = 0.621371

def GCDistance(lat1,lon1,lat2,lon2):
	lat1*= DE2RA
	lat2*= DE2RA
	lon1*=DE2RA
	lon2*=DE2RA
	d = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1 - lon2)
	return AVG_ERAD * acos(d)

def GCAzimuth(lat1, lon1, lat2, lon2):
	result = 0.0
	ilat1 = int(0.50 + lat1 * 360000.0)
	ilat2 = int(0.50 + lat2 * 360000.0)
	ilon1 = int(0.50 + lon1 * 360000.0)
	ilon2 = int(0.50 + lon2 * 360000.0)
	lat1*= DE2RA
	lat2*= DE2RA
	lon1*=DE2RA
	lon2*=DE2RA
	if ilat1 == ilat2 and ilon1 == ilon2: 
		return result
	elif ilat1==ilat2:
		if ilon1<ilon2:
			result = 90.0
		else:
			result = 270.0
	elif ilon1==ilon2:
		if ilat1 > ilat2:
			result = 180.0
	else:
		c = acos(sin(lat2)*sin(lat1) + cos(lat2)*cos(lat1)*cos((lon2-lon1)))
		if c==0: c=1
		A = asin(cos(lat2)*sin((lon2-lon1))/sin(c))
		result = (A * RA2DE)
		if ilat2<ilat1 and ilon2!=ilon1:
			result=180.0-result
		elif ilat2>ilat1 and ilon2<ilon1:
			result+=360.0
	return result

def ApproxDistance(lat1, lon1, lat2, lon2):
	lat1*= DE2RA
	lat2*= DE2RA
	lon1= -DE2RA * lon1
	lon2= -DE2RA *lon2
	F = (lat1 + lat2) / 2.0
	G = (lat1 - lat2) / 2.0
	L = (lon1 - lon2) / 2.0
	sing = sin(G)
	cosl = cos(L)
	cosf = cos(F)
	sinl = sin(L)
	sinf = sin(F)
	cosg = cos(G)

	S = sing*sing*cosl*cosl + cosf*cosf*sinl*sinl
	C = cosg*cosg*cosl*cosl + sinf*sinf*sinl*sinl
	W = atan2(sqrt(S),sqrt(C))
	R = sqrt((S*C))/W
	H1 = (3 * R - 1.0) / (2.0 * C)
	H2 = (3 * R + 1.0) / (2.0 * S)
	D = 2 * W * ERAD
	return D * (1 + FLATTENING * H1 * sinf*sinf*cosg*cosg -FLATTENING*H2*cosf*cosf*sing*sing)

def EllipsoidDistance(lat1, lon1, lat2, lon2):
	distance = 0.0
	faz = 0.0
	baz = 0.0
	r = 1.0 - FLATTENING
	if lon1==lon2 and lat1==lat2: return distance
	lat1*= DE2RA
	lat2*= DE2RA
	lon1*=DE2RA
	lon2*=DE2RA

	cosy1 = cos(lat1)
	cosy2 = cos(lat2)

	if cosy1 == 0.0: cosy1 = 0.0000000001
	if cosy2 == 0.0: cosy2 = 0.0000000001
	c=0.0
	d=0.0
	e=0.0
	cz=0.0
	c2a=0.0
	sa=0.0
	sx=0.0
	cx=0.0
	y=0.0
	sy=0.0
	cy=0.0

	tu1 = r * sin(lat1) / cosy1
	tu2 = r * sin(lat2) / cosy2
	cu1 = 1.0 / sqrt(tu1 * tu1 + 1.0)
	su1 = cu1 * tu1
	cu2 = 1.0 / sqrt(tu2 * tu2 + 1.0)
	x = lon2 - lon1

	distance = cu1 * cu2
	baz = distance * tu2
	faz = baz * tu1
	while(fabs(d - x) > EPS):
		sx = sin(x)
		cx = cos(x)
		tu1 = cu2 * sx
		tu2 = baz - su1 * cu2 * cx
		sy = sqrt(tu1 * tu1 + tu2 * tu2)
		cy = distance * cx + faz
		y = atan2(sy, cy)
		sa = distance * sx / sy
		c2a = -sa * sa + 1.0
		cz = faz + faz
		if c2a > 0.0: cz = -cz / c2a + cy
		e = cz * cz * 2. - 1.0
		c = ((-3.0 * c2a + 4.0) * FLATTENING + 4.0) * c2a * FLATTENING / 16.0
		d = x
		x = ((e * cy * c + cz) * sy * c + y) * sa
		x = (1.0 - c) * x * FLATTENING + lon2 - lon1
	x = sqrt((1.0 / r / r - 1.0) * c2a + 1.0) + 1.0
	x = (x - 2.0) / x
	c = 1.0 - x
	c = (x * x / 4.0 + 1.0) / c
	d = (0.375 * x * x - 1.0) * x
	x = e * cy
	distance = 1.0 - e - e
	distance = ((((sy * sy * 4.0 - 3.0) * distance * cz * d / 6.0 - x) * d / 4.0 + cz) * sy * d + y) * c * ERAD * r
	return distance

def getKmPerLonAtLat(dLatitude):
	dLatitude *= DE2RA
	return 111.321 * cos(dLatitude)

def getLonPerKmAtLat(dLatitude):
	return 1 / getKmPerLonAtLat(dLatitude)

def getKmPerLat():
	return 111.000

def getLatPerKm():
	return 1 / getKmPerLat()\

def ConvKilometersToMiles(dValue):
	return dValue / 1.609344

def DMStoDEC(deg,min,sec):
	"""Converts DMS ( Degrees / minutes / seconds ) to decimal format longitude / latitude"""
	return deg+(((min*60)+(sec))/3600)

def DECtoDMS(dec):
	"""Converts decimal longitude / latitude to DMS ( Degrees / minutes / seconds )"""
	deg=int(dec)
	tempma = (dec - deg) * 3600
	min = floor(tempma / 60)
	sec = tempma - (min*60)
	return (deg,min,sec)